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DEFINITION OF SYMBOLS 


Symbol 

A 

n 

A 2, 2 

B 2 

c 

n, m 
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g e 

GM 
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n, m 


n, n 


m 


P (sin d>) 
n 

P (sin <6) 
n 


Definition 

coefficients in the series expansion for the geoid radius 

coefficient of the sectorial harmonic in the geoid expansion 

coefficient in the series expansion of the geoid radius 

numerical coefficients in the geopotential function 

flattening or oblateness 

ellipticity (or flattening) of the equator 

gravity at the equator 

mean equatorial gravity 

gravitational parameter of the earth ( = m) 

coefficients of zonal harmonics in geopotential 

coefficients of tesseral harmonics in geopotential 

coefficients of sectorial harmonics in geopotential 

oblateness coefficients ( older notation) 

centrifugal factor at the equator 
also: index in series notation 

associated Legendre polynomials of degree n and order m 
Legendre polynomials of nth degree 
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Definition 

equatorial radius of the earth 
mean equatorial radius of the earth 
polar radius of the earth 
north polar radius of the earth 
south polar radius of the earth 
longest equatorial radius of the earth 

shortest equatorial radius of the earth 
radial distance from the center of the earth 
potential function 

geographic longitude, positive east of Greenwich 

longitude of the longest equatorial radius (semi-major 
axis) 

gravitational parameter of the earth 

factor in equatorial gravity acceleration equation 

geocentric latitude 

factor defined by equation ( 22) 

angular velocity of the earth's rotation 

centrifugal factor at the equator 

perturbation parameter in potential function 
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SURVEY AND COMPARATIVE ANA LYS IS OF 
CURRENT GEOPHYS I CAL MODELS 

SUMMARY 


By using currently accepted constants of geodesy, especially several 
different sets of oblateness coefficients, seven models of the geometry of the 
earth and its gravitational field are numerically derived and compared with each 
other and with a newly developed geopotential based on recent measurements of 
the secular and periodic perturbations of artificial satellite orbits. 

Conclusions are drawn which show, among other things, that the currently 
used NASA Standard Model deviates significantly from the "best" geoid, and that 
there are better analytical, as well as numerical, models by which improved 
approximation can be achieved. 

Also, it is shown that improvements in the numerical values of the used 
oblateness coefficients, rather than in their number, appear to have little influ- 
ence on the values of the gravity, but that they affect the geometries significantly. 
Increasing the number of zonal harmonics coefficients, however, can lead to 
considerable improvements in the gravity model, at least in the range studied 
( up to and including J^) . 


INTRODUCTION 


In the past a number of theories of varying complexity have been developed 
to describe the figure and the gravitation field of the earth for use in geodetic, 
astrodynamic, and astronomical calculations. 

In scrutinizing these theories, a number of discrepancies are manifest 
which appear to be caused by two factors; i.e. , ( 1) inherent inconsistencies in 
the development of some of the theories, and (2) differences in gravitational 
and geometrical parameters used by their authors. 

With regard to the first consideration, while the earth's gravitational 
potential (i.e. , the potential energy of the earth in relation to the position 
relative to the earth) is commonly expressed by the potential function as a 



series of spherical harmonics, as adopted by the International Astronomical 
Union ( IAU) [ 1 ] , the shape of the earth introduced in the geopotential at the 
earth's surface with the simultaneous assumption that the potential along the sur- 
face of the geoid is constant, is that of an ellipsoid which approximates the actual 
geoid. The flattening of the ellipsoid is usually determined approximatively by 
the astrogeodetic methods of geometrical geodesy (sometimes also by the less 
accurate techniques of gravimetry) , using variants of Clairaut's theorem. How- 
ever, even if higher-order extensions of the latter are used, they remain only 
approximations, which become increasingly inaccurate for increasing deviation 
of the geopotential from that of a regular ellipsoid. 


The second consideration refers to the fact that, with our increasing 
knowledge of our globe and its gravity anomalies, the descriptive geographical 
parameters are subject to change, particularly the number and values of the 
oblateness coefficients, but also — to a lesser degree — the gravitational param- 
eter GM q , and the geometrical parameters flattening f, and mean equatorial 

radius R e . Urgently needed standardization of the principal parameters was 

accomplished by NASA in 1963, but new and additional data have been determined 
continuously ever since and are being determined from an ever-increasing popu- 
lation of artificial satellites in precision-tracked ( Baker-Nunn) orbits of all 
inclinations . 


To relate these new observations to the geoid without having to accept 
the inconsistencies of the second- and third -order Clairaut extensions, a new 
theory was developed recently at MSFC by H. Krause [2,3] , which further in- 
creased the bulk of existing geographic parameters and models. 

It is the objective of the following analysis to define and to compare 
seven different geographic models, one of which is the current (but outdated) 
NASA Standard. These models are based on four different sets of oblateness 
coefficients, all of which are applicable to present practical work, and also on 
the fact that for each set of oblateness conditions two body geometries can be 
defined; i.e. , a spheroid composed of superimposed spherical harmonics and a 
regular ellipsoid approximating this spheroid. 

In the first part of the analysis, the general background and the nature 
of the problem is considered in some detail, followed by the derivation and 
description of the geophysical theories. In the third part, the models are 
compared with each other. A number of salient conclusions, forming the fourth 
part, can then be draw. 
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GENERAL BACKGROUND AND DISCUSSION OF THE PROBLEM 
The Potential Function of the Earth 


In the study of geocentric motion such as that encountered on rocket 
trajectories and satellite orbits, among the most important constants necessary 
for describing the environment are those pertaining to the gravity field acting on 
the moving body, and to the shape or figure of the gravitating body, the earth. 

Theories attempting to describe the earth's gravity field are, by necessi- 
ty, fundamentally based on Newton's Universal Law of Gravitation, which is 
rigorously accurate for a central force field as produced by a perfectly homo- 
geneous and perfectly spherical body. The gravitational potential of such a 
Newtonian force field is 


U = - =■ 
r 


( 1 ) 


The potential function of the earth U^, or the geopotential, is defined 

as the integral of the gravity forces over the entire field, which are the resultant 
of the gravitational forces due to Newtonian attraction and the centrifugal force 
due to the rotation of the earth and its atmosphere. Thus, the geopotential is 
the sum of the potential of the gravitational field and the potential of the centrif- 
ugal force. 


U @ = -^-(l + iM+ -|ca 2 r 2 cos 2 <£ (2) 

If the rotational term is omitted, equation ( 2) expresses the gravitational 
potential which is appropriate to inertial coordinates, as useful for exo-atmos- 
pheric (free-flight) flight phases. For geodesy, as well as for all other applica- 
tions requiring a rotating earth coordinate system, the term is retained. 

In equation ( 2) , the modifier ip accounts for the fact that, contrary to 
the basic assumptions of the Newtonian potential, the earth's internal mass dis- 
tribution is not entirely homogeneous and its shape is not spherical, but oblate. 
Also, its surface is covered with irregularities, such as continental highlands, 
depressions, mountain ranges, ocean deeps, valleys, etc., which give rise to 
gravitational anomalies and require an additive perturbative term to complete 
the geopotential . 
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Since the potential function of the earth is a solution of the Laplace 
equation (V 2 U^ = 0) , it is a harmonic function and can be represented ana- 
lytically as the sum of a series of two-dimensional (spherical) harmonics in a 
manner similar to the one -dimensional Fourier series as representation of a 
nonanalytical function on a circle. The coefficients of the terms of this infinite 
(converging) series are polynomials of the general Legendre type, while the 
terms themselves are the so-called surface harmonics, the introduction of 
which is due to Legendre and Laplace. If the spherical harmonic is a function 
of two variables, such as latitude and longitude (rather than only one, the lati- 
tude) , it involves the so-called Associated Legendre Polynomials. 

The acceleration of gravity g of the earth is defined as V U and can 
therefore also be expressed in spherical harmonics. 

The perturbative parameter ip of the gravitational potential in equation 
(2) can be expanded, as described, in an infinite series of spherical harmonics. 


°o n / R \ n 

ib = / / 1 ] P (sin <f>) [C cosmA + S sinmA 1 , 

^ ^ r / n v ^ n, m n,m 1 

n=l m=0 (3) 


where r is the distance from the center of the earth, R is the earth’s mean 

e 

equatorial radius, <±> is the latitude, A is the longitude, C and S are 
M 5 ^ n, m n,m 

m 

numerical coefficients, and the P^ are the associated Legendre polynomials. 

Introducing equation (3) in equation (2) , one obtains the general formula for 
the earth potential, as recommended by the IAU (1962) [4, 1], 


U = — 


- 

oo 

n / 

R \ 



v 1 

e 

1 + 

L 

n=l 

3 \ 

m=0 

77 




P m (sin cb) ( C cos mA 
n n, m 


+ S sinmA ) 
n, m 


1 9 9 9 

+ — co r cos 0 


(4) 


The Spherical Harmonics 


It is seen that the second (aspherical) part of equation (4) consists of 
an infinite series of harmonics which — in a rigorous sense — describes the 
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perturbative potential accurately only if an infinite number of terms are con- 
sidered. Taking the terms by families according to the indices, it is further- 
more seen that for the first existing family, where m = 0, i < n < °o , the equa- 
tion reduces to 


U 


m=0 


a 

r 



P (sine*) 
n 



(5) 


for the purely gravitational potential. 

In this special case, P^(sin <p) are the (simple) Legendre polynomials 

of argument sirup , usually referred to as M zonal" harmonics. Obviously, the 
potential U ^ is independent of the longitude A; hence, it is constant along 

(latitude) parallels. This means that the undulations of the harmonics on the 
sphere are zero along parallels and become alternately positive and negative 
between these n parallels in a distance ir . Thus, equation (5) refers to an 
axially symmetric earth which is fluted along the parallels and where the sym- 
metry axis is the rotational axis through the poles. 


use 


As an alternate for C 

n, 0 


it has been recommended by the IAU [1] to 


J 

n 



(6) 


Thus, for the body-of-revolution case, the purely gravitational potential 
becomes 


- 

CO 

/ R 

n 

i 

1 

- X J ' 1 

n=l n 

e ' 

1 • P (sirup) 


(7) 


The true figure of the earth does not exactly exhibit rotational symmetry, 
as assumed by setting m=0 in equation (4) . The principal axes of inertia of the 
true earth do not exactly coincide with the rotational axis system. Thus, the 
angular momentum of the earth is not pointing exactly in the same direction as 
the angular velocity. Consequently, products of inertia do exist, which are 
represented by higher harmonics. By letting m assume integral values from 
1 to n, terms are added to the second part of equation (5) , which include the 
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longitude A as a second variable. Consequently, the associated Legendre poly- 
nomials P m (sind>) are also required, as well as coefficients J . These 
n n, m 

terms are referred to as tesseral harmonics , with the exception of the two 
terms in equation (4) resulting from letting m = n, so that 

C P n cos nA and S P n sinnA 
n, n n n, n n 

appear. These two terms are called "sectorial" harmonics . 

The tesseral harmonics ( 0 < m < n) have zeros both along meridians 
and along parallels; thus, the undulations represented by them become zero on 
the reference spheroid simultaneously on a number of meridians and parallels. 
They undulate in the form of a network across the spheroid similar to a chess- 
board composed of fields that are alternatingly positive and negative. Thus, 
for example, the spherical harmonic represented by the function Pg (with the 
coefficient J 6> 4 *) forms a network of five meridians (four zeros in 180 degrees) 
and (6-4) = 2 parallels (not equator) on the sphere, cutting it into 30 (alter- 
natingly raised and depressed) pieces. 

In the special case of the two families of sectorial harmonics, the undu- 
lations of the function become zero on the reference spheroid for cos <fi = i 1 
(on the poles) and for cos n A = 0, i. e. , along n meridians which are symmetri- 
cal to n planes through the rotational axis, thus cutting the spheroid in sectors 
which are alternately concave and convex. J 4 4 , for example, would represent 
undulations featuring four meridians and eight sectors. 


1. The notation used commonly in the literature for the tesseral and sectorial 

J-coefficients is J . Thus, the first sectorial harmonics coefficient would be 
nm 

J 2 2 , the first tesseral coefficient J 12 , etc. While this practice represented an 

ambiguity in notation, with respect to the zonal harmonics, no great confusion 

could be caused as long as the zonal harmonics used in practical work did not 

exceed the tenth order. With present analyses of the higher zonal harmonics 

already going to n = 21, m = 0 [5] however, notation of the mentioned type would 

easily confuse tesseral with zonal harmonics. Thus, the zonal J 14 could easily 

be interpreted as the tesseral n = 1, m = 4. Since powers of J are sometimes 

n 

encountered, the old practice of writing is also not desirable. Consequently, 

it is suggested that tesseral and sectorial coefficients are indexed with a dividing 

comma, i.e., J , while the zonal notation J be retained. The notation 
n, m n 

P^ for the associated Legendre polynomials should also remain unchanged. 
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The Harmonics Coefficients 


The coefficients of the tesseral harmonics of the earth, especially the 
higher ones, are not really known. The sectorial coefficients J 2 , 2 and J 4) 4 
have been determined tentatively during the past 6 years and are sometimes 
used in high precision orbit determination programs [6] . For J 2> 2 , which has 
two zeros along a parallel (four sectors) and thus determines the ellipticity 
of the equator and the parallels, Kaula [7] gives 

J 2j2 = (1.80 ± 0 . 1) x 10~ 6 . 

Of the tesseral harmonic perturbations of satellite orbits, only the 
sectorial J 2 2 is large enough to be of practical concern in most orbit analyses; 
specifically, for example, this harmonic affects the supplemental energy require- 
ments of 24-hour satellites. The sectorial harmonic due to J 2? 2 is depicted in 
Figure 1. 

Because of the uncertainty in and relative insignificance of the tesseral 
harmonics coefficients, they are usually neglected in the potential function which 
therefore assumes the commonly used form 



r 00 

/ R \ n 

_ 

U = £ 

1 - ^ J • 


■ P (sin<i) 

r 

L n=l n 

v r ) 

n 


which was first suggested by Brouwer [8] and adopted by the IAU in 1961. If 
so desired, tesseral and sectorial terms with the appropriate associated 
Legendre polynomials can be added to the second part of the equation, as required 
by the inclusion of the longitude A in the analysis, as well as the rotational 
potential ( 1/2) r 2 c o 2 cos 2 <p to the whole equation to transform the gravitational 
potential in inertial coordinates for a free body to the gravity potential for a 
body attached to the earth. 

Of the ( theoretically infinitely many) zonal harmonics coefficients J^, 
only a few are known with adequate accuracy. If an equatorial geocentric 
coordinate system is assumed as a reference system such that its origin coin- 
cides with the center of mass, and if the latitude 0 is measured relative to 
this center, the first term in the sum of equation ( 8) , which is proportional to 
1/r , vanishes since the coefficient J l5 multiplied with R^, would represent the 

distance between coordinate origin and center of mass; thus, Jj= 0 . 
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( meters ) 



FIGURE 1. FIRST SECTORIAL HARMONIC (DUE TO J 2> 2 ) AT EQUATOR 

Of the other zonal functions, the second order harmonic is by far the 
largest, since it refers to the mass distribution due to the earth oblateness, 
thus defining the departure of the shape of an earth ellipsoid from that of 
an ideal sphere, while J 4 being three orders of magnitude smaller than J 2 , and 
all following even J's account for deviations in hemispherical mass distribution 
symmetrical about the equator plane, i.e. , for deviations of the geoid from the 
ellipsoid. Each even J of increasing order, while maintaining symmetry about 
the equator plane, "narrows down" the figure into an increasingly complex 
shape, while the sum of subsequent even J's describe the deviation of the geoid 
from the preceding shape. The even J’s, however, all represent spheroid 
shapes with identical northern and southern "hemispheres." 
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The asymmetry of the northern and southern "hemispheres" of the geoid, 
also referred to as "pear" shape, is taken into account by the odd surface har- 
monics, which contain odd powers of siruj > , especially by the principal coefficient, 
J 3 . Higher-odd-order J's are at least one order of magnitude smaller than this 
coefficient. 

The first zonal harmonic is shown in Figure 2 2 . Figures 3 and 4 depict 
the zonal harmonics due to J 3 and J 4 , and J 5 and J 6 , respectively. Of particular 
interest is the ovoid form of J 3 , which is reponsible for the "pear" shape of the 
earth . 


The Theorem of Clairaut 

In the past, higher order zonal harmonics, and especially the odd har- 
monics, have been neglected in the geopotential. This was mainly due to the 
methods of determining the coefficients used before the advent of the artificial 
earth satellites. Also, since the amplitudes of high harmonics — according to 
the quantity 1/r in their terms — decay very rapidly with altitude above the sur- 
face, and their effect on bodies in the earth's external gravitation field dwindles 
rapidly with increasing order and/or altitude, coefficients above fourth order, 
as well as the odd term J 3 , were usually neglected. Before the advent of the 
artificial satellites, the moon was the only moving body in the earth's potential, 
but at its distance only the second harmonic is still discernible in its motion. 

Determinations of the J-coefficients J 2 and J 4 were fundamentally based 
on geodetic surveys, gravity measurements, observations of the moon, and 
Clairaut's theorem. The latter deals with the form of a surface which encloses 
all the matter of a rotating body with various density distributions, which con- 
stitutes an equipotential surface (surface of constant potential) . Clairaut's 
formula (1743) is 

J = f - 2 m ’ < 9 ) 

where J is the second harmonics coefficient, m is the centrifugal factor at the 
equator w 2 R/g (withcu 2 R= centrifugal acceleration at equator, g = gravita- 

0 O 

tional acceleration at equator) , and f is the flattening of the ellipsoid of revolu- 
tion, as assumed here. 


2. The figure is actually the second term of equation (34) . 


9 






North Pole 



South Pole 


FIGURE 3. THIRD AND FOURTH ORDER SPHERICAL HARMONICS 
TERMS [16 P 3 (sirup ); 20 P 4 (sin0)] IN EARTH GEOID 
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North Pole 



South Pole 


FIGURE 4. FIFTH AND SIXTH ORDER SPHERICAL HARMONICS 
TERMS [ 1 . 3 P 5 (sin0) ; -4 P 6 (sin0) ] IN EARTH GEOID 
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Obviously, Clairaut’s formula relates the amplitude of the second har- 
monic in gravity to that of the corresponding harmonic in the radius vector of 
the equipotential surface. 


Clairaut’s formula, because of the approximations made in deriving it, 
is not accurate enough for practical use, A second order extension of the theorem 
has been derived by Darwin ( 1899) , DeSitter ( 1924) , Jeffreys ( 1954) , Bullen 
( 1946) , Cook ( 1959) , et al. , resulting in equations involving higher orders of 
the flattening, such as the third-order equation given by Kaula [1], 


J 2 




|m- |f+ fm 2 + i|f 2 + 0(f 3 ) 


( 10 ) 


With the use of Clairaut's theorem, as well as its extensions by the 
mentioned authors, the earth flattening — before the artificial satellites — could 
be derived separately in two ways: from a harmonic analysis of observed gra- 
vity values (where the measured gravity anomalies are used to compute the 
absolute deflections of the vertical and, thus, the true undulations of the geo- 
potential surface; also, to derive the detailed gravity field) and from astrono- 
mical observations of the constant of precession and the theory of the moon's 
motion (DeSitter, 1915) , leading to J 2 . Higher J's could then be computed from 
the higher-order extensions of Clairaut's theorem, using the derived flattening, 
but when the first artificial satellites were flown, it was found that higher-order 
harmonic coefficients determined with this method did not agree with the J's 
actually observed from the orbital perturbations of the satellites. The reason 
for this serious discrepancy is simply the fact that for the higher-order terms, 
due to the irregular density distribution within the earth, the basic assumption 
of an ellipsoid of revolution as the figure of the earth in the Clairaut formula 
and all higher-order extensions is an increasingly bad approximation of the geoid. 


The Geoid 

Distinction is made between the geodetic geoid, which depends on the 
surface irregularities of the earth, and the geophysical geoid, which is governed 
by the internal irregular mass distribution of the earth. 

The geodetic geoid is the true figure of the earth. While the actual earth 
surface with its irregularities is not an equipotential surface, the geophysical 
geoid, however, is required to be an equipotential level, identical with mean 
sea level and its continuation under the continents. The geoid, thus, is an 
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extremely complex shape, requiring a very large number of higher-order har- 
monics for its description. 

The reason for the use of the geoid in geodesy is the following: In 
principle, the geometrical form of the earth's surface can be found independent 
of observations of the potential and gravity, by geodetic triangulation methods, 
and then, given the values of these quantities on the surface, Laplace's equation 
can be solved for all the exterior space. Because of the extreme complexity 
of the true boundary, this would be an almost hopeless task. It is therefore 
universal practice to refer all geodetic observations to the gravitational equi- 
potential surfaces and to determine the form of these from geodetic observations. 
Because the sea-level equipotential surface is an internal surface on land, the 
gravity field is not computed for the actual earth but for a model earth which is 
related to the actual earth and which is bounded by an equipotential surface. 

The field so computed will not agree with that of the actual earth throughout 
space, but the model may be chosen so that it agrees where observations can be 
made. This is the geoid. 

The geodetic geoid must be computed point by point and cannot be given 
by few parameters, since it is not a mathematical surface but depends on the 
irregular distribution of visible and invisible masses near the earth's surface. 
For geodetic surveys, however, this point-by-point calculation would be imprac- 
tical. Therefore, reference spheroids of revolution have been used instead. 

The reference spheroid is assumed to be an equipotential surface of the same 
volume and flattening as the geoid. 

The history of geodesy has seen the following principal reference 
spheroids (all ellipsoids) : 

Germany (and several European states) . . . .1841 . . .Bessel 


England 1880. . .Clarke 

U. S.S.R 1938. . .Krassowski 

U.S 1866 . . .Clarke 

U.S 1910 . . .Hayford 

U.S 1963. . . Kaula/Fischer 



Determination of the Harmonics Coefficients 


As has been mentioned, J 2 in former times was determined from the 
moon's motion. For the lower amplitudes of the higher order harmonics, how- 
ever, these effects are not discernible because r becomes too large. However, 
when the first artificial satellites were orbited, it was found that the gravita- 
tional anomalies expressed by the higher-order harmonics characterize the 
potential field at their altitudes. Thus, the existence and magnitude of J 3 
(i.e. , the "pear" shape of the earth) was first discovered from satellite 1958/32 
(Vanguard-1), launched in March 1958 [9]. 

The asphericity of the earth, represented by the spherical harmonics, 
causes constant secular and periodic perturbations on a satellite orbit. By 
observing the secular and periodic perturbations, the zonal harmonics coeffi- 
cients can be determined. Since the secular changes, especially the orbital 
precession (i.e., the regressive rotation of the nodal line) and the advance of 
the perigee (i.e. , the rotation of the apsidal line) , depend primarily on the 

even-numbered J , while the long-period oscillations in four of the five orbital 
n 

elements (eccentricity, inclination, longitude of ascending node ( right ascension), 
and argument of perigee) are caused primarily by the odd harmonics [10], the 
determination of the even is usually done from secular changes, while the 

odd are computed from periodic perturbations. An orbital perturbation 

suffered by one particular satellite yields one linear equation between the 
harmonics. For example, if /3 is the amplitude of the observed oscillation in 
eccentricity, we obtain an equation for the coefficients of the odd harmonics, 

J 3 , J 5 , J 7 . . . of the form [5] 


AsJ 3 + a 5 j 5 + A 7 J 7 + . . . - j3 , (ii) 

where the A's are constants for a particular orbit and depend mainly on the 
orbital inclination i and semi-major axis a. By using several different satel- 
lites, several equations of the form (11) can be obtained. Therefore, if a 
large number of satellites (widely distributed over as many inclinations as 
possible) is available for analysis, an equally large number of can be deter- 
mined. 


In solving the simultaneous equations ( 11) , the main error of the analy- 
sis is introduced by the number of to be determined remaining finite and all 

terms higher than the highest desired J being assumed negligible. In fact, 
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however, these higher coefficients are not necessarily negligible, and since 
their effect is — by necessity — lumped together in the computed J^, in any 

one analysis of a given set of satellites there are always several different solu- 
tions to the J's, depending on the number of J's considered in the algebraic 
equations. Thus, for example, King-Hele, et al. [5], recently computed odd 
zonal harmonics up to J 21 from eccentricity oscillations of 14 satellites, giving 
two different sets, i.e. , one with J 3 to J 15 , and the other with J 3 to J 21 . 


Also, the equations are sometimes solved by several differend methods 
to obtain more confidence in the results. For example, the equations may be 
solved for three J's at a time, assuming the higher J’s to be zero, in fours and 
fives, and also with least-squares and minimum -residual methods. The 
coefficients given by King-Hele [5] in the larger set are 


J 3 = 

(-2. 50 ± 0. 01) 

X 

1 

0 

'I — 1 

J 5 = 

(-0. 26 ± 0. 01) 

X 

10~ 6 

J 7 = 

(-0. 40 ± 0. 02) 

X 

10“ 6 

J 9 = 

(0 ± 0. 06) x : 

10“ 6 


Jil = 

(-0.27 ± 0.06) 

X 

10“ 6 

Jl3 = 

(+0.36 ± 0. 08) 

X 

10“ 6 

J« = 

(-0. 65 ± 0. 10) 

X 

10~ 6 

Jj7 = 

(+0.30 ± 0. 08) 

X 

10“ 6 

J 19 “ 

(0 ± 0. 11) x 10“ 6 


J21 = 

(+0.58 ± 0. 11) 

X 

10“ 6 


( 12 ) 


I 

The latest most complete set of zonal harmonics, both even and odd 
coefficients, has been evaluated by Kozai. In 1963, Kozai [11] derived a set 
of values for eight coefficients of zonal harmonics (up to and including J 9 ) from 
the then available observations of secular motions of the node and the perigee 
and the amplitudes of long-periodic terms of artificial satellites. Later, in 
1964, from precisely reduced Baker-Nunn observations of nine high-inclination 
satellites (28 degrees through 95 degrees) he produced [12] a new set of 
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thirteen coefficients (up to J 14 ) : 

J 2 = (1082.645 ± 0.006) X 10 “ 6 


J 3 

= 

(-2.546 ± 

0 . 020 ) 

X 

10 " 

J 4 

= 

(-1. 649 ± 

0.016) 

X 

10 “ 

J 5 

= 

(- 0.210 ± 

0.025) 

X 

10 " 

J 6 

= 

(+0. 646 ± 

0 . 030) 

X 

10 “ 

J 7 

= 

(-0. 3 33 ± 

0. 039) 

X 

10 " 

J 8 

= 

(-0.270 ± 

0. 050) 

X 

10 “ 

J 9 

= 

(-0. 053 ± 

0.060) 

X 

10 " 

J 10 

= 

(-0. 054 ± 

0. 050) 

X 

10 “ 

Jll 

= 

(+0. 302 ± 

0.035) 

X 

10 " 

Jl2 


(-0.357 ± 

0. 047) 

X 

10 “ 

J 13 

- 

(-0. 114 ± 

0. 084) 

X 

10 “ 

J 14 

= 

( + 0. 179 ± 

0. 063) 

X 

10 “ 



(13) 


The values of Kozai's even coefficients have been compared recently [13] 
with two other sets of King-Hele/Cook and Smith. The agreement is excellent 
for all inclinations greater than 28 degrees, which was the smallest inclination 
of available satellite orbits. It is therefore concluded that, as long as there 
are no satellites at inclinations between 10 degrees and 25 degrees the evalua- 
tion of the even harmonics in the potential carried beyond the present status 
does not seem worthwhile. For the odd harmonics, this is not the case, since 
their primary effect, the long-periodic changes (e.g., the amplitude of the 
oscillation in eccentricity) decrease to zero as the inclination tends to zero. 

For this reason, King-Hele felt justified to publish the odd zonal harmonics to 
J 21 , as given above, with an indication that J 23 and J 25 are small. 
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Krause's Theory of the Geoid Surface 


If Kozai's value for J 2 , 


J 2 = (1082.645 ± 0.006) X 10" 6 


is introduced in a relation of a third -order Clairaut theory [1], 


* - - h 


3 

2 m 



lif + 

49 



(14) 


a value of f = 1/298. 254 results for the flattening of the reference earth ellip- 
soid. Using this f to compute the theoretical value of J 4 for the reference ellip- 
soid assumed to be in hydrostatic equilibrium, one obtains -2.350 x 10~ 6 , as 
compared to the observed value of -1.649x 10 -6 . This discrepancy illustrates 
the increasing inaccuracy of Clairaut's (ellipsoid) theorem and its extension, 
if applied to higher -order harmonics. 


To find a more accurate relation between the radius vector of the earth 
and its gravity field, thereby cleaning up the inconsistency in the older theories, 
Krause recently developed a new theory [2,3] which essentially introduces the 
radius vector of the geoid, rather than of an ellipsoid, consistently expressed 
to the same accuracy as the potential field it is related to. Since the new theory 
is kept general, it can handle all higher order zonal harmonics up to with the 

same internal consistency, which the older theories, based on Clairaut, Darwin, 
Helmert, et al. , are unable to do. Using Krause’s theory with Kozai's thirteen 
coefficients and introducing the radius vector of the geoid expanded in spherical 
harmonics up to J 14 in the similarly expanded geopotential, a new flattening of 
the earth can be computed 3 , f = 1/298. 1840, which differs both from the old 
third -order value given above and from the flattening of the Kaula/Fischer 
ellipsoid used currently as a world datum, for example, in trajectory computa- 
tions ( 1/298.30) . It is believed that the value derived with Krause's theory is 
more consistent with the observed spherical harmonics coefficients. 

In his theory [2,3] , Krause has made use of the fact that the generalized 
( equipotential) surface of the geoid, similarly to the geopotential, can be repre- 
sented alternatively by two series expansions, one consisting of spherical har- 
monics, the other of powers of the sine of the latitude. When introduced in the 
geopotential, the chosen series must be of the same degree as the potential 


3. The value given in references 2 and 3 has been slightly corrected by the 
present writer. 
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function used. The author shows that the coefficients of the terms of the series 
(A n in the spherical harmonics series, in the power series) are related to 

the known oblateness coefficients J of the geopotential equation U and can be 

n 

computed from these. Insertion of the radius equation thus obtained in the 
equation of the total gravity acceleration in the direction of the normal yields 
the gravitational acceleration of the geoid as an expansion of surface harmonics, 
based on a consistent geoid geometry. Because of the relation between the 

( or B n ) coefficients and the oblateness coefficients J^, and assuming the surface 

to be of constant potential, the oblateness of the northern and southern "hemis- 
pheres" and, thus, the mean meridional flattening f can be computed from the 

A coefficients, relating now to the true geoid and the full set of observed J , 
n n 

rather than to a reference ellipsoid, as in Clairaut's theory. 

For the geoid radius, expanded in spherical harmonics, 


OO 

-R 

— — = A Q - 2j A P (sin<£) + 3A 2) 2 cos 2 <p cos 2( A - X fl ) , (15) 

R n=l n n 

e 

in which the first sectorial harmonic, with the associated Legendre polynomial 
P 2 (sin0) exactly given by 3 cos 2 cp , is taken into account to express the ellip- 
ticity of the earth equator, Krause finds for the A -coefficients [2,3] : 


OO 



(16) 

(17) 

(18) 
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(k* 0,2,4) 


(19) 


A 

K 


- J 

X K 


where 






GM 


centrifugal factor at equator, 


(20) 


B, = - 


r v 1 • 3 • 5. . . ( 2v + 3) 1 ~ 

v J 2„ + 2 + 2“e 

L „,o 2 v\ 2'. 


E<-‘> 


( 21 ) 


and 


X = 1 - 


E <-‘> 

n=l 


n 3 . 5 • 7. . . (2n + 1) 
2 • 4 • 6. . . 2n 


n ~ CO 

2n e 


( 22 ) 


If the first variation of U with longitude, as represented especially by 
the first sectorial term, is taken into account as it was in equation ( 15) , where 
3 cos ^cp is used as an exact substitute of P| (sin<£) , the geopotential of the earth 
becomes in inertial coordinates 


U 


GM 

r 


E ■ 

n=2 


n 



P ( sirup ) 
n 


+ 



cos 2 cb cos 2( A - A ) 
o 


+ 


~ 2 j 

— CO cos ^ 
2 e ^ 


(23) 


Since the three components (in polar cooi’dinates ) of the gravity are 
defined as 


8 U 

^r 8 r 

1 8 U 

r cos p 8 A 


(24) 
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( 25 ) 


i 9 u 

^ <p r 3 <p 

while the total gravity acceleration in the direction of the normal is 

g m ^ 4 + 4 + ’ 


the gravity, including the sectorial term, can be given by determining the poten- 
tial derivatives and expanding the sum of their squares according to equation 
(25) in a series approximating the square-root of equation (25) . The result is 
[2,3] 


g = 



oo 


2 <■> + 

n=2 


1) J 


n 



P^sin <p ) 


+ 9 J2,2 



cos 2 <b cos 2( A - A ) 

'r o 






3 

COS 2 4 ) 


+ 



+ 



sin 2 cp cos 2 cp 


i 

j 


(26) 


where only the derivative of P 2 ( sin (p) with respect to <p has been included 
[expressed as 

P 2 (sirup) = — = 9 sirup cos cp ] , 

while the higher order derivatives have been neglected. 


ANALYSIS OF MODELS 
Model 1 (J2 _ 0nly) 

For computation of trajectories close to earth, featuring small central 
range angles, the effects of the equatorial bulge and other irregularities can 
usually be considered negligible, primarily because the flight times involved 
are relatively short. In this model, it is therefore assumed that the earth is 
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an ellipsoid of revolution, symmetrical around all three axes. Except for the 
second harmonic, all higher order effects of asphericity (which are three or 
four orders of magnitude smaller than J 2 ) are taken as zero. 


With these assumptions, the gravity equation (26) reduces to 


g 




P 2 (sin4>) 



+ 



. 9 

sin 


<P COS 2 (t> 


\ 

\ 


(27) 


It would be more accurate to compute the gravity directly from the 
potential of the ellipsoid, rather than from the series approximation, equation 
( 27) . The task of determining the potential of a homogeneous ellipsoid was a 
very formidable one and remained unsolved till the end of the 18th century, 
despite attempts at a solution by Newton, Maclaurin, D'Alembert and Lagrange. 
Finally, Laplace succeeded in finding the general solution, followed later by 
Gauss and others. The general form of the potential of a homogeneous, triaxial 
ellipsoid for an exterior point involves an elliptic integral of the first kind; 
hence, the components of the gravity force of the ellipsoid also contain elliptic 
integrals. For the present case of the ellipsoid of revolution, the equations 
become considerably simpler; the elliptic integrals reduce to logarithmic - 
cyclometric integrals. The investigations of the above mentioned authors 
generally refer to homogeneous bodies. Clairaut was the first who investigated 
inhomogeneous bodies. 


(28) 


( 29 ) 


A third-order relation of Clairaut's theory is given by Kaula [ 1 J : 


= l f l 1 - i f 


i 

3 m 


1 - ^m - |f+ T ni 2 + 777 f 2 + 0 ( f 3 ) 
2 74 49 v ’ 


Kaula gives also 


GM = 3 . 9 8 6 0 32 x 10 20 cm 3 sec -2 

© 


398603. 2 km 3 sec -2 


co = 0. 720, 211, 585 x 10 4 sec 1 


R = 6378165.0 m 
e 


g^ = 978. 0300 cm -sec' - 


-2 
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Using these values, the centrifugal factor at the equator, m, becomes 


co 2 • R 

m = = 0.003,467,773,255, 

g e 


(30) 


which is in agreement with reference 7. 

For Model 1, the value for J 2 used is the first one of the set given 
recently by Kozai [ 12 ], equation ( 13) : 

J 2 = 1082.645 x 10" 6 . (31) 

In a strict sense, this J 2 -value is not completely consistent with the 
above given value of g^ from Kaula [1] . The equatorial gravity of the oblate 

spheroid may be computed from 

gm © 

B e = a • It • 

e 


where the factor cr accounts for two effects. First, it includes a correction to 
Mq because the mass of the earth’s atmosphere, included inM^, does not 

contribute to the gravitational acceleration at the surface, and second, it modi- 
fies the succeeding term, valid for a sphere only, for the oblateness effects. 

In the Kaula value of equation (29) , a by necessity refers to the older J 2 and 
probably also to the older J 4 . However, the inconsistency is not felt to be of 
practical consequence for the present purpose. For the J 2 of equation (31) , 
a would become 0.9981616959 (by solving equation (22) and computing a ~ X ( 1”A), 

A = M /M ~ = 0. 8594 x 10“ 16 [3]) , so that the above equation results in 

atm © 

g^ = 978. 026 gal, 4 mgal lower than Kaula's value, which is presently accepted 
as standard [ 14] . 


It must be kept in mind that the assumption of all higher being zero 

introduces, in the strict sense, an error in the model which is not taken into 
account by the above J 2 -coefficient, since its particular value has been deter- 
mined from simultaneous algebraic equations involving 12 additional coefficients 
of higher order, which were different from zero. However, it is assumed that 
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the uncertainty caused by this fact will be of no concern for the present pur- 
poses. This assumption is justified by the results of this analysis. 

Introducing the J 2 -value in equation (28) , one obtains the oblateness 

f = 0.0033528465 = 1 : 298.254 , 
which is consistent with the assumed geopotential. 

With the value for R given in equation (29) , the polar earth radius 
becomes 

R = R (1 - f) (32) 

p e ' 

R = 6356780.0 m 
P 

R = 6378165.0 m 
e 

and the general radius at latitude 0 

R R 

R = ' - - 6 P (33) 

N R cos 2 0 + R sin 2 0 
p ^ e 

Since the assumption of an ellipsoid is not fully consistent with the 
gravitation assumed (based on J 2 ) , a second body has been determined by using 
the first spherical harmonic in a truncated expression of the type ( 15) . By 
solving equation ( 16) , ( 17) , ( 20) , (21) and (22) , the coefficients are deter- 
mined: 

A 0 = 0. 9988810301 
A t = 0 

A 2 =- 0.2236813952 x 10~ 2 

A 3 = 0 

A 4 = - 0. 1501014034 x 10 -5 
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Neglecting the coefficient A 4 , the radius of the spheroid is then deter- 
mined from 

Rj = 6371028.026 - 14266.76846 • P 2 (sin<j p) , (34) 

where P 2 (sin<£) is the second-degree Legendre polynomial. 

The resulting spheroid is depicted in Figure 5. 

The total gravity, equation (27) , and the surface shape of the ellipsoidal 
Model 1, equation (33) , as well as the surface of the body formed by the J 2 
harmonic, equation (34), are given in Table I. 


Model 2 (NASA Standard) 

For trajectories of longer duration, such as those of upper stages, as 
well as for near-earth orbits, the geopotential exhibits anomalies such that its 
effect deviates from that of an attracting ellipsoid of revolution. For almost 
all practical purposes, however, it suffices to add only the two next higher 
order zonal harmonics to the potential and to consider the gravity anomalies 
averaged over all longitudes, thus retaining the rotational symmetry around the 
polar axis and rendering tesseral and sectorial terms unnecessary. 

The second model is represented by the geopotential and the astro- 
geodetic world datum as adopted by the Ad Hoc NASA Standard Constants 
Committee [14] in 1963 for Project Apollo and other NASA programs. The 
standard is based on the Kaula/Fischer Ellipsoid of 1963. Irene Fischer of 
the U. S. Army Map Service, in 1960, published a world ellipsoid based on an 
imposed flattening f = 1/2 98. 30 (which at that time appeared to be the best 
available value from early artificial satellites) , by determining the ellipsoid 
of revolution which best fit the geoid contours derived from astrogeodetic 
measurements [15]. Kaula later [1], in 1963, published a new value for the 
equatorial radius, 

R e - 6378165. 0 i 25. 0 m , (35) 

which was a compromise between Fischer’s value [15] and a 1961 value by 
Kaula [ 16] . The NASA Committee meeting at Goddard Space Flight Center in 
May 1963 adopted this value, along with the following constants [14] : 
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TABLE I. MODEL i 


Latitude 

Radius 

Radius 

Height of 

Gravity 

North or South 

( Ellipsoid) 

( Spheroid) 

Spheroid over 

(gal) 

(deg) 

(m) 

(m) 

Ellipsoid (m) 


90 (poles) 

6356780. 0 

6356761.3 

-18.7 

983.2049 

85 

6356941. 6 

6356923. 8 

1 

<1 

00 

983. 1663 

80' 

6357421. 7 

6357406. 6 

-15. 1 

983. 0516 

75 

6358205. 8 

6358194. 8 

-11.0 

982. 8641 

70 

6359270. 5 

6359264.6 

-5.9 

982. 6093 

65 

6360583. 7 

6360583. 5. 

- 0.2 

982.2944 

60 

6362106. 1 

6362111.3 

5. 2 

981. 9286 

55 

6363791. 7 

6363801.7 

10. 0 

981. 5228 

50 

6365589. 7 

6365603. 3 

13.6 

981.0890 

45 

6367445. 6 

6367461.3 

15. 7 

980.6401 

40 

6369303. 1 

6369319. 4 

16. 3 

980. 1897 

35 

6371105. 7 

6371121.0 

15. 3 

979. 7517 

30 

6372798. 5 

6372811. 4 

12. 9 

979. 3394 

25 

6374329. 7 

6374339.2 

9. 5 

978. 9656 

20 

6375652. 3 

6375658. 1 

5. 8 

978. 6422 

15 

6376725. 7 

6376727. 9 

2.2 

978.3794 

10 

6377517. 0 

6377516. 1 

- 0.9 

978. 1854 

5 

6378001. 7 

6377998. 9 

- 2.8 

978. 0664 

0 (equator) 

6378165. 0 

6378161. 4 

- 3.6 

978. 0300 


GM _ = 398603.2 km 3 sec -2 

© 

f = 1 : 298.30 
g = 978. 030 gal 

G 

J = 1. 62345 (±0. 00030) x 10 


-3 


H = -0. 575 (±0. 025) x 10' 


-5 


D = 0. 7875 (±0. 0875) x 10 

R = 6356783.3 m 

P 

co = 3461. 414, x 10" 6 sec -1 

e a 


(36) 
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where the J, H, D coefficients apply to an alternate (older) form of the geopo- 
tential, and the value 

w 2 R 3 
~ e 

W e " GM 

is fully consistent with the assumed parameters. 

In the modern notation, the oblateness coefficients become 

j 2 = f- J = 1082.30 x 10 -6 | 

i 3 I 

j 3 = |-H = -2.30X10" 6 ) (37) 

D ( 

J 4 = - = - 1. 80 x 10“ 6 I 

It is pointed out that the above set of constants is consistent within 
itself for most practical purposes. However, since the geometrical figure 
assumed is still an ellipsoid of revolution (equal northern and southern 
"hemispheres") , derived from J 2 through a third-order relation of Clairaut's 
theory, equation (28) [ 1 ], it is at best only an approximation and contains 
theoretically an inconsistency. 

Table II lists the gravity of Model 2 (second column) , computed from 
equation (26) without the longitude -dependent term, with J 2 , J 3 , J 4 as given by 
equation (37) , and in the third column the local radius of the associated Kaula/ 
Fischer reference ellipsoid, computed from equation (33) with the data given 
in (35) and (36) . 

As an illustration of the inherent inconsistency of the NASA Standard 
model, it is seen that the gravity field at the North Pole is higher by 18.4 mgal 
than at the South Pole, while the radial distance computed from the model is 
unable to account for this "pear" shape. 

To compare the deviation of the ellipsoid from the (here assumed) geoid, 
the shape of the geoid must be computed. The new theory by Krause [ 2 , 3] , 
provides a relatively simple method for doing this. To write the geoid radius in 
an expansion of the second, third and fourth surface harmonics, the coefficients 
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TABLE II. MODEL 2 



Latitude 

Gravity 

Radius 

Radius 

Difference 


(degrees) 

(gal) 

Ellipsoid 

Geoid 

(m) 




(m) 

(m) 



90N (pole) 

983.2231 

(1356733. 3 

6356796.3 

-13.0 


85 

983. 1839 

6356944. 9 

6356957.7 

-12. 8 


80 

983.0676 

6357424. 9 

6357437. 1 

-12.2 


75 

982. 8776 

6358208. 9 

6358220. 1 

-11.2 


70 

982.6195 

6359273. 4 

6359283. 2 

- 9.8 

N 

65 

982.3011 

6360586.4 

6360594. 6 

- 8.2 

O 

60 

981. 9317 

6362108. 6 

6362114. 9 

- 6.3 

R 

55 

981.5225 

6363793. 9 

6363798.3 

- 4.4 


50 

981. 0857 

6365591 . 6 

6365593. 9 

- 2.3 

T 

45 

980. 6346 

6367447.2 

6367447. 6 

- 0.4 

H 

40 

980. 1 *29 

6369304. 5 

6369303. 1 

1.4 


35 

979. 7444 

637 1 106. s 

6371104. 0 

2. 8 


30 

979. 3325 

6372799. 4 

6372795. 4 

4.0 


25 

97s. 959s 

6374330 . 3 

6374325. 7 

4.6 


20 

97s. fi3 Ml 

63 75ti 52. « 

6375648 . 0 

4. 7 


13 

97 s . 3770 

63 76726 . 0 

637672 1 . 7 

4.3 


10 

97s. 1 s,5u 

6 3 7 i . 1 1 i . 1 

6377513. 9 

3.2 


5 

97 s . 06? s 

637M101 . s 

637 8000 . 0 

1.8 

0 (equator) 

97."> . 0.300 

637M 65. 0 

637 si 65. 0 

0 


5 

978. 0701 

6378001. 8 

(>378003. 8 

- 2.0 


10 

978. 1895 

6.3775 1 7. 1 

6377521. 1 

- 4.0 


15 

978.3833 

6376726. 0 

6376731. 9 

- 5.9 


20 

978. 6455 

6375652. 7 

6375660. 1 

- 7.4 


25 

978. 9679 

6374330. 3 

6374338. 8 

- 8.5 

S 

30 

979.3404 

6372799.4 

6372808.3 

- 8.9 

O 

35 

979.7514 

6371106. 8 

6371115. 4 

- 8.6 

U 

40 

980. 1884 

6369304. 5 

6369311. 9 

- 7.4 

45 

980. 6379 

6367447. 2 

6367452. 8 

- 5.6 

T 

50 

981 . 0862 

6365591.6 

6365594. 7 

- 3. 1 

H 

55 

981. 5198 

6363793. 9 

6363794.0 

- 0. 1 

60 

981 . 9257 

6362108. 6 

6362105. 4 

3.2 


65 

982.2919 

6360586. 4 

6360579. 9 

6. 5 


70 

982. 6074 

6359273. 4 

6359263.7 

9. 7 


75 

982. 8629 

6358208. 9 

6358196. 4 

12. 5 


80 

983. 0509 

6357424. 9 

6357410.3 

14.6 


85 

983. 1660 

6356944. 9 

6356929. 0 

15. 9 


90S (pole) 

983. 2047 

6 3 56 7 S3. 3 

6356766 . 9 

16.4 
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ill I 

in equation (15) must be determined with J 2 , J3, J 4 , equation (37) , from 
equations (16) through (22). 

One obtains 

A 0 = 0.9988805375 

Aj = 0 

A 2 = 0.2236440156 x 10" 2 * * 
A 3 = - 0.2304227311 x 10 -5 
A 4 = - 0. 3313107217 x 10" 5 

with which the geoid surface of Model 2 can be expressed by 

R 2 = 6371024.88 - 14264.38 • P 2 (sin0) + 14.69 • P 3 (sin0) 

+ 21.13 • P 4 (sin0) (m) (39) 

where the P n (sin<£) are the Legendre polynomials of argument sin<p . 

Column 4 of Table II lists the results of equation (39) , while the devia- 
tions of the ellipsoid from the geoid are given in column 5. It is seen that the 
NASA Standard Model results in an earth radius which is 13m below the geoid at 
the North Pole, 16. 4m above the geoid at the South Pole, and between zero and 
12m above or below it in the intermediate latitudes. 



With 


f 


R pN + R pS 


2 R 


the mean meridional oblateness becomes 


f = 0.0033525944 = 1 : 298.276, 

differing slightly from the Standard value of 1 : 298.30 . 


(40) 


(41) 
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The NASA Standard, here called Model 2, adopted by the Ad Hoc Com- 
mittee in 1963, has served and is serving its purpose well, which is to ensure 
accuracy and consistency between the various NASA agencies and contractors 
participating in the space program. "While during 1961-63 a real attempt was 
made to select a set of constants which might be termed the "best available at 
the time" [14], the chief qualification of the adopted set was standardization. 
Since the adoption of the NASA Standard for Apollo and other space programs, 
new and better determinations of the geographical (and astrodynamic) constants 
have been made. The Standard of 1963, still valid today, is therefore not based 
on more current information; however, with the present status of the Apollo 
program, it would probably be ill-advised to update the Standard at this time. 

The following models make use of the more recent determinations in 
various combinations. 


Model 3 (New J3, J4) 

Similar to the preceding analysis, a model is assumed which takes the 
"pear" shape of the earth into account as well as the second even zonal harmonic, 
but not any dependence of the geopotential on longitude. The spherical harmonics 
coefficients are taken from Kozai [12], equation (13), 

J 2 = 1082.645 x 10 " 6 

J 3 = - 2.546 x i(T 6 , (42) 

J 4 = - 1.649 x 10 -e 


where they have actually been determined together with ten additional coefficients, 
which are here assumed zero. This introduces a slight error of no practical 
concern for the present purpose. 

Gravity is again computed from equation (26), without the sectorial 
term, with J 2 , J 3 , J 4 as given in (42) . 

According to the new J-coefficients, the a -factor in the equation for the 
equatorial gravity acceleration must now be slightly adjusted; g^ becomes now 

g = 978. 029 gal 
e 

instead of the Standard value 978. 030 gal. 
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The reference ellipsoid (of revolution) postulated for this model is 
required to be consistent with the geopotential used. 

Using Krause's theory, the A-coefficients, equation (16) through (22), 
in equation (15) are obtained. 

A 0 = 0. 9988804206 
A t = 0 

A 2 = 0.2236788408 x 10" 2 

A s = - 0.2550678850 x 10~ 5 

A 4 = - 0.3160571410 x 1(T 5 
and the radius of the assumed geoid becomes 

R 3 = 6371024.14 - 14266. 61 P 2 (sin0) + 16. 268 P 3 (sin^) 

+ 20. 16 P 4 (sin0) (m) . (44) 

The results of this equation for the geoid surface are listed in column 4 
of Table III. It can be seen that the polar radii are 

R = 6356794. 0 m 
pN 

R _ = 6356761.4 m 
pS 

so that the mean meridional flattening, equation ( 40) , becomes 

f = 0.003353206 = 1 : 298.222 . 

Comparison with the geoid of Model 2 shows that the new values for 
J 2 , J 3 , J 4 cause the Model 3 geoid to be smaller at most latitudes (by 2.3 m at 
the North Pole, 5. 5 m at the South Pole, and by 1.2 m and 0. 6 m at 45N and 
45S, respectively) . 

For the reference ellipsoid, the equator radius imposed is the currently 
accepted value due to Kaula, equation (35) . With the above flattening, the polar 
radius of the ellipsoid then becomes 
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TABLE m. MODEL 3 


Latitude 

(deg) 

Gravity 

(gal) 

Radius 

Ellipsoid 

(m) 

Radius 

Geoid 

(m) 

Difference 
(Height of 
Ellipsoid 
(m) 

90 

983.2239 

6356777.7 

6356794.0 

-16.3 

85 

983.1848 

6356939.4 

6356955.4 

-16.0 

80 

983.0685 

6357419.5 

6357434. 8 

-15.3 

75 

982.8786 

6358203. 7 

6358217. 9 

-14.2 

70 

N 

^ 65 

982.6205 

6359268.4 

6359281.2 

-12.8 

982.3020 

6360581.8 

6360592.8 

-11.0 

° 60 

981.9326 

6362104.3 

6362113.3 

- 9.0 

R 55 

981.5234 

6363790.1 

6363796. 8 

- 6.7 

T 50 

981.0866 

6365588.3 

6365592.6 

- 4.3 

T 45 

980.6354 

6367444.4 

6367446.4 

- 2.0 

H 

40 

980. 1835 

6369302. 1 

6369302.0 

0. 1 

35 

979.7449 

6371105.0 

6371103.0' 

2.0 

30 

979.3328 

6372797.9 

6372794. 5 

3.4 

25 

978.9601 

6374329.3 

6374324. 9 

4.4 

20 

978. 6382 

6375652.0 

6375647.3 

4.7 

15 

978.3771 

6376725.6 

6376721.2 

4.4 

10 

978. 1851 

6377516.9 

6377513.5 

3.4 

5 

978.0679 

6378001.7 

6377999. 8 

1.9 

0 

978.0294 

6378165.0 

6378165.0 

0 

5 

978.0705 

6378001.7 

6378004.0 

- 2.3 

10 

978. 1900 

6377516. 9 

6377521.5 

- 4.6 

15 

978.3840 

6376725.6 

6376732.4 

- 6.8 

20 

978.6464 

6375652.0 

6375660. 8 

- 8.8 

25 

978. 9689 

6374329.3 

6374339.4 

-10. 1 

30 

c 

979.3416 

! 6372797.9 

6372808. 8 

-10. 9 

o 

o 35 

979.7527 

6371105.0 

6371115.6 

-10.6 

° 40 

980.1896 

6369302.1 

6369311.8 

- 9.7 

U 45 

980.6389 

6367444. 4 

6367452.2 

- 7.8 

T 50 

981.0871 

6365588.3 

6365593.4 

- 5.1 

T 55 
H 

981.5205 

6363790. 1 

6363792. 1 

- 2.0 

60 

981.9261 

6362104.3 

6362102.7 

1.6 

65 

982.2919 

6360581.8 

6360576.5 

5.3 

70 

982.6070 

! 6359268.4 

6359259.6 

8.8 

75 

982.8623 

6358203.7 

6358191.8 

11.9 

80 

983.0501 

6357419.5 

6357405.2 

14.3 

85 

983. 1650 

6356939.4 

6356923.6 

15.8 

90 

983.2037 

6356777.7 

6356761.4 

16.3 
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R = R (i-f) = 6356777.70 m , 
p e 

and the ellipsoid shape can then be computed from equation (33) . Results are 
listed in column 3 of Table in. The deviations of the ellipsoid from the geoid 
are given in column 5. The largest excursions are again at the poles, with 16.3 
m height difference. 


Higher-Order Bodieswith Equator Symmetry 


Model 4a (Triaxial El lipsoid) . In the preceding three ellipsoids, the 
longitude dependence of the geopotential has not been taken into account. For 
the present model, it is assumed that the geoid can be better approximated by 
a triaxial ellipsoid. How well it agrees with the geoid is one of the objectives 
of this investigation. 

In a triaxial ellipsoid, the equator and all parallel sections exhibit a 
certain ellipticity, expressed — in the geopotential — by the sectorial coefficient 
J 2; 2 * The body is symmetrical with respect to the equator plane; hence, the 
odd-order terms in the spherical harmonics series of the geopotential (more 
exactly, of the radius vector of the surface) are zero. Thus, the earth is now 
assumed to be no longer symmetrical around the rotational axis. 


The general formula of a triaxial ellipsoid is 

r.2 


_x 

R 


2 yi_ 

+ j>2 


R 2 


= i 


'0 


-l 


(45) 


where R and R are the two equatorial semi-axes. By introducing polar 
e 0 ej 

coordinates, Krause [3] has derived a general formula for the radius vector, 


R 

e 

R 


cos *<$ 


1 + 7 f2 
4 e 


f cos 2 (X - A ) 
e o 



sin 2 <ft 

( 1 - f) 


% 


(46) 


where R g is the mean equatorial radius, determined by the astrogeodetic and 
gravimetric methods of geodesy [17, 18] and currently assumed to be [i] 
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I 


R = 6378165. 0 m . 
e 

The symbol 0 is the geocentric latitude of the parallel, A is the longitude of 
the point on the surface, measured positive east of Greenwich, f is the merid- 
ional flattening of the ellipsoid, f is the ellipticity of the equator (and all 

0 

parallel sections) , and A q is the longitude of the apsidal line of the equatorial 

ellipse (i.e. , the angle between the radius vector to the reference meridian on 
the equator and the semi-major axis of the equator ellipse, positive east of 
Greenwich) , as well as the longitude of the associated parallel ellipses in both 
"hemispheres." Kaula [7] gives A q = - 14.5 degrees ±1.5 degrees. In Ref- 
erence 1, he quotes -8 degrees to -25 degrees; however, more recent data 
indicate an improved value of [3] 

A q = - 18 degrees ± 3 degrees (47) 

Thus, the semi-major axis of the equator ellipse is rotated by 18 degrees 
west of the Greenwich meridian. Consequently, the largest equator radius is 
located approximately 600 statute miles north of Ascension Island in the Mid- 
Atlantic, and — diametrically opposed — in the vicinity of the Solomon Islands 
in the Pacific, while the shortest radii are approximately 500 statute miles 
south and west of Ceylon in the Indian Ocean and about 1800 statute miles west 
of Peru in the Pacific Ocean. 


By solving Krause's radius equation, equation (15), twice for A = A q 

and A = A ±90 degrees, respectively, at the equator (0=0) and introducing 
o 

the derived expressions for R and R in the definition of the flattening (or 
oblateness) , 0 1 


R 


R 


(48) 


R 


one obtains a simple expression for the ellipticity of the equator, viz. , 


f = 
e 


6 * A 2j 2 = 



(49) 
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The best current value for J 2> 2 is given by Kaula [7] : 

J 2> 2 = (i. 8± 0. 1) x 10" 6 . (50) 

Using the most recent set of oblateness coefficients, equation (13) , in 
the equation for the factor x> equation ( 22) , one obtains 

X = 0.9981691858 . (51) 

With x and J 2j 2 , equation (50), equation (49) yields 

f = 10.81980906 x 10 -6 = 1:92423 . (52) 

e 

From equation ( 48) it is then found for the difference between the semi- 
major and semi -minor axis of the equator 

AR = R - R = 69. 0 m . (53) 

e e 0 ej 

With the mean equatorial radius, R , according to equation (35), the 
longest radius of the equator becomes 

_ AR 

R + — - 6378199.5 m 

e 2 

and the shortest radius 
_ AR 

R - — = 6378130. 5 m 

e 2 

at the geographical locations described above. 

The sectorial harmonic is depicted in Figure 1. 

The flattening f of the triaxial ellipsoid must be the same as for the 
"best" geoid (Model 5a and 5b) treated later, since the flattening equation, 
equation (58), contains only even zonal harmonics, in accordance with the 
symmetry behavior of the latter. With the even J of equation ( 13) , equation 
( 58) yields 
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f = 0.003353633015 = 1 : 298. 184 . 


The polar radii then become 

R = (1 -f)R = 6378775.0 m . 
p e 

The surface of the triaxial ellipsoid has been computed for various 
latitudes and longitudes from equation (46) . Results are given in Table IV, 
column 2. 

The gravity of the Model 4a ellipsoid can be determined from equation 
( 26) . All available even are used, equation ( 13) , with the odd zonal coeffi- 
cients set equal to zero according to the assumed symmetry with respect to 
the equator plane. With x as given in equation (51) , the factor a becomes 
0. 9981683280, and the mean gravitation at the equator is then 


g 


e 




978. 0320 gal 


differing by 2 mgal from Kaula's Standard value. 

Results are presented in Table IV, column 5. 

Mo del 4b (Symmetrical Geoid) . With the oblateness coefficients up to 
J 14 readily available, it is of interest to investigate the deviation of the above 
triaxial ellipsoid from the surface of a quasi -ellipsoid, or rather "symmetrical 
geoid, " derived from a series expansion of the known even zonal harmonics 
equivalent to the above employed gravity expansion. Using Krause's theory and 
solving equations (16) through (22), the coefficients inequation (15) are found 
to be 


A 0 = +0.9988800610 
A 2 = + 0.2236730320. x -10 -2 
A 4 = -0.3180618021 x 10” 5 
A 6 = - 0.6471848752 x 10 -6 
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TABLE IV. MODELS 4a AND 4b 



Longitude 

Radius of 

Radius of 

Height of 

Gravity of 

Latitudes 

( deg) 3 

Tri -ellipsoid 

Symm. Spheroid 

Spheroid 

Spheroid 


EOG 

WOG 

(m) 

(m) 

(m) 

(gal) 

90N/90S 



6356775. 0 

6356775.0 

0 

983. 1815 

(Poles) 







75N/75S 

0 

180 

6358203.0 

6358203.0 

0 

982.8385 


30 

150 

6358200.9 

6358200. 9 

0 

982.8376 


60 

120 

6358199.0 

6358199.0 

0 

982.8367 


90 

90 

6358199.3 

6358199.3 

0 

982.8368 


120 

60 

6358201.4 

6358201.4 

0 

982.8378 


150 

30 

6358203.2 

6358203.2 

0 

982.8386 

60N/60S 

0 

180 

6362109.2 

6362114.2 

4.9 

981.9077 


30 

150 

6362101.4 

6362106.3 

. 4. 9 

1 981.9041 


60 

120 

6362094. 5 

6362099.3 

4.8 

981.9009 


90 

90 

6362095.4 

6362100.2 

4.8 

981.9013 


120 

60 

6362103.2 

6362108. 1 

4.9 

981.9049 


150 

30 

6362110.1 

6362115.1 

4. 9 

981.9081 

45N/45S 

0 

180 

6367456.9 

6367461.7 

4. 7 

980.6259 


30 

150 

6367441.2 

6367446.0 

4.7 

980.6186 


60 

120 

6367427. 4 

6367432.0 

4.6 

980.6122 


90 

90 

6367429.2 

6367433. 8 

4.6 

980.6131 


120 

60 

6367444.8 

6367449. 6 

4.7 

980.6203 


150 

30 

6367458.7 

6367463.5 

4.8 

980.6267 

30N/30S 

0 

180 

6372818.1 

6372819.3 

1. 1 

979.3374 


30 

150 

6372794. 6 

6372795.6 

1.0 

979.3263 


60 

120 

6372773.7 

6372774. 7 

1.0 

979.3166 


90 

90 

6372776.4 

6372777.4 

1.0 

979.3179 


120 

60 

6372800.0 

6372801.0 

1.0 

979.3288 


150 

30 

6372820. 8 

6372822.0 

1.2 

979.3384 

15N/15S 

0 

180 

6376751.4 

6376849 .3 

-2. 1 

978.3883 


30 

150 

6376722.0 

6376719.9 

-2.1 

978.3747 


60 

120 

6376696.0 

6376693.9 

-2. 1 

978.3627 


90 

90 

6376699.4 

6376697.2 

-2.2 

978.3643 


120 

60 

6376728.7 

6376726.6 

-2. 1 

978.3778 


150 

30 

6376754. 8 

6376752.7 

-2.1 

978.3898 

Equator 

0 

180 

6378192. 9 

6378192. 9 

0 

978.0448 


30 

150 

6378161.4 

6378161.4 

0 

978.0303 


60 

120 

6378133. 5 

6378133.5 

0 

978.0175 


90 

90 

6378137. 1 

6378137.1 

0 

978.0192 


120 

60 

6378168.6 

6378168. 6 

0 

978.0337 


150 

30 

6378196.5 

6378196.5 

0 

978.0465 


a. EOG = East of Greenwich 
WOG = West of Greenwich 



Ag = - 0.2704952265 X 

O 

i 

05 

(54) 

A 10 = - 0. 5409904530 x 

h-^ 

O 

1 


A 12 = - 0.3576547995 x 

10~ 6 


A 14 = + 0. 1793283168 x 

10~ 6 



The surface equation, equation (15), then becomes 

R 4 = 6371021.84 - 14266.235 • P 2 ( sin0) + 20. 29 • P 4 (sin 0) 

-4.13 • P 6 (sin0) + 1.73 • P 8 (sin 0) + 0.345 • P 10 (sin<£) 

+ 2.28 • P 12 (sin0) - 1.14 • P 14 (sin0) 

+ 34.51 • cos 2 0 • cos 2( X + 18°) (m) . (55) 

Equation (55) has been computed for some values of <fi and X ; results 
are given in column 3 of Table IV. The deviation of the symmetrical geoid from 
the true triaxial ellipsoid is listed in the fourth column of Table IV. Agreement 
is seen to be extremely good. Height differences are around 5m in the higher 
latitudes. The largest height difference occurs around 55 degrees North or 
South latitude, with over 6m for all longitudes. 

"Best" Geoid Bodies 

Model 5a ("Best" Geoid with Sectorial Term) . The most realistic model 
of the earth and its potential field must obviously make use of all available 
spherical harmonics coefficients, obtained from real-world observations. 

Using the theory developed recently by Krause [2, 3], it is — as shown above — 
relatively simple to express the geoid surface, assumed an equipotential level, 
in an expansion of spherical harmonics and appropriate coefficients. This "true” 
geoid is assumed as Model 5a of this investigation. 

For the basic imposed constants (mean earth radius at equator, rate 
of the earth’s rotation with respect to inertial space, and the gravitational 
parameter of the earth) the best currently available values [1] are chosen: 

R = 6 378 165. 0 m 
e 

co = 0.729211585 x 10 -4 sec -1 
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398603.2 ± 3. 0 km 3 sec -2 


GM © " 
and for the additional parameters 

J 2}2 = ( 1-8 ± 0 . 1 ) x 10- 6 

A = - 18° ± 3° . 
o 


Consistent with the above values, the mean equatorial gravitation is 
§ = 978.0320 gal (see Model 4a) 

Solving again equations (16) through (22) for all available J including 

odd orders, equation (13), the following A -coefficients of equation (15) are 
obtained: n 


A o 

= 

+ 0. 9988800610 



Ai 

= 

0 



A 2 

= 

+ 0. 2236730329 

X 

10 -2 

A 3 

= 

- 0.2550669802 

X 

10 -5 

A 4 

= 

- 0.3180618021 

X 

10" 5 

A 5 

= 

- 0. 2103851762 

X 

10 -6 

A 6 

= 

+ 0.6471848752 

X 

10-® 

A 7 

= 

- 0.3336107793 

X 

10“ 6 

A 8 

= 

- 0.2704952265 

X 

10- 6 

A 9 

= 

- 0. 5309721113 

X 

10“ 7 

A 10 

= 

- 0. 5409904530 

X 

10- 7 

A 11 

— 

+ 0.3025539200 

X 

10 -6 




(56) 
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A 12 = - 0.3576547995 

X 

10 -6 

Aj3 = - 0.1142090956 

X 

10 _e 

A 14 = + 0.1793283168 

X 

io - 6 


with 

X = 0.9981691858 

B 2 = 0.3383437858 x 10~ 2 


/ 


The equation for the surface of the geoid becomes then 4 

R 5 = 6371021. 844 - 14266.235 P 2 (sin0) + 16.268 P 3 (sin0) 

+ 20.286 P 4 (si n<p) + 1. 34 P 5 (sin0) - 4. 127 P 6 (si n<£) 

+ 2.127 P 7 (sin<£) + 1. 725 P 8 (sin0) + 0. 338 P 9 (sin<£>) 

+ 0.345 P lo (sin0) - 1.929 P^sin cp) + 2.28 P 12 (sin0) 

+ 0. 728 P 13 (sin$) - 1 . 144 P 14 ( sin<£) 

+ 34. 505 cos 2 <p cos 2(A + 18°) , (m) (57) 

where P^(sin<£) are the Legendre polynomials of argument sin<fi 5 . 

By writing the radius equation, equation ( 15) , for the poles ( <fi = ± 90 
degrees) , one equation each results for the northern and the southern "hemi- 
spheres, M which since 


f 


R + R 
pN pS 


2 R 


lead to an expression for the mean meridional oblateness of the geoid [2, 3] 


4. Some of the amplitudes (coefficients) of the geoid-equation in reference 
3 have been found to be slightly incorrect. The values are corrected here. 

5. See Appendix A. 
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.V=i 


i - (-D 


v 1 • 3 • 5. . . (2u - 1) 
2 • 4 • 6. . . 2v 


2v 


+ K + B *(jh - " e ' B *)} 

where x anb B 2 are given in equation ( 56) and co = 3461. 414 x 10 -6 sec -1 . 

G 


( 58 ) 


With Kozai's harmonics coefficients up to J 14 , equation (58) yields 
f = 0.003353633015 = 1 : 298. 184, 

which is different from the currently accepted (Standard) value of 1960, viz., 


f = 1 : 298.30 , 

which is consistent with the old J 2 , equation (37) , over a third-order Clairaut- 
type expression only, as mentioned earlier. 

The surface of the "best" geoid, equation (57), and its gravity, equa- 
tion ( 26) , have been computed and are tabulated in Table V for various latitudes 
and longitudes. It can be seen that the geoid, in accordance with the fact that 
of all tesseral harmonics only the first sectorial was included, is symmetrical 
around a plane defined by the major axis (at A = - 18 degrees) of the equator 

ellipse and the axis of rotation. If higher tesserals were taken into account, 
this last symmetry would also disappear. 

The sectorial harmonic at the equator is shown in Figure 1. The curve 
depicts (greatly exaggerated) the height of the "best" geoid over the rotationally 
symmetric geoid of Model 4b around the equator. At latitudes 45N and 45S, 
the amplitudes of the harmonic, for equal meridians, are half their value at 
the equator. 

Because of the sectorial term depicted, the longest equator radius 
(also the longest radius vector to all latitudes) is at 162 degrees east and at 
18 degrees west of Greenwich, with 6 378 199. 5 m. The shortest radial vector 
is at 72 degrees east and 108 degrees west of the prime meridian, with 
6 378 130. 5 m. The difference amount to 69. 0 m. 
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TABLE V 


GEO ID 


Latitudes 

a 

Longitude 

Radius of 
Geoid 

Gravity of 
Geoid 

Longitude a 

Radius of 
Geoid 

Gravity of 
Geoid 


EOG 

WOG 

(m) 

(gal) 

EOG 

WOG 

(m) 

(gal) 

90N 

(North 

Pole) 



6356793.9 

983. 18807 





75N 

0 

180 

6358217.6 

982. 84373 

90 

90 

6358213. 8 

982.84200 


30 

150 

6358215.4 

982. 84275 

120 

60 

6358215.9 

982. 84298 


60 

120 

6358213.6 

982.84189 

150 

30 

6358217.8 

982. 84384 

6 ON 

0 

180 

6362118.1 

981.90820 

90 

90 

6362104.2 

981.90178 


30 

150 

6362110.3 

981. 90457 

120 

60 

6362112.1 

981.90540 


60 

120 

6362103.3 

981.90136 

150 

30 

6362119.0 

981.90861 

45N 

0 

180 

6367458.7 

980.62507 

90 

90 

6367430. 8 

980.61223 


30 

150 

6367443.0 

980.61782 

120 

60 

6367446.6 

980.61948 


60 

120 

6367429.0 

980.61140 

150 

30 

6367460.5 

980.62590 

3 ON 

0 

180 

6372812.7 

979.33540 

90 

90 

6372770. 8 

979.31614 


30 

150 

6373789. 0 

979.32452 

120 

60 

6372794.4 

979.32701 


60 

120 

6372768. 1 

979.31489 

150 

30 

6372815. 4 

979.33664 

15N 

0 

180 

6376743.6 

978.38634 

90 

90 

6376691.6 

978.36238 


30 

150 

6376714.2 

978.37281 

120 

60 

6376721.0 

978.37591 


60 

120 

6376688.2 

978.36083 

150 

30 

6376747.0 

978.38789 

0 i 

0 

180 

6378192. 9 

978.04484 

90 

90 

6378137. 1 

978.01916 

( Equator) 

30 

150 

6378161. 4 

978.03034 

120 

60 

' 6378168.6 

978.03366 


60 

120 

6378133.5 

978.01750 

150 

30 

6378196.5 

978.04650 

15S 

0 

180 

6376755.0 ! 

978.39020 

90 

90 

6376702. 9 

978. 36623 


30 

150 

6376725.6 

978.37667 

120 

60 

6376732.3 

978.37976 


60 

120 

6376699.6 

978.36469 

150 

30 

6376758.4 

978.39174 

3 OS 

0 

180 

6372825. 9 

979.33889 

90 

90 

6372784.0 

979.31963 


30 

150 

6372802.2 

979.32801 

120 

60 

6372807.7 

979.33050 


60 

120 

6372781.3 

979.31838 

150 

30 

6372828.6 

979.34013 

45S 

0 

180 

6367467. 7 

980.62674 

90 

90 

6367436. 8 

980.61390 


30 

150 

6367448. 9 

980.61949 

120 

60 

6367452. 5 

980.62115 


60 

120 

6367435.0 

980.61307 

150 

30 

6367466. 5 

980.62757 

60S 

0 

180 

6362110.2 

981. 90722 

90 

90 

6362096.3 

981. 90080 


30 

150 

6362102.3 

981. 90360 

120 

60 

6362104. 1 

981. 90443 


60 

120 

6362095. 4 

981. 90039 

150 

30 

6362111. 1 

981. 90764 

75S 

0 

180 

6358188.5 

982.83336 

90 

90 

6358184.7 

982. 83164 


30 

150 

6358186.3 

982.83239 

120 

60 

6358186.8 s 

982. 83261 


60 

120 

6358184.5 

982. 83153 

150 

30 

6358188.7 

982. 83347 

90S 

(South 

Pole) 



6356756. 1 

983. 17496 






a. EOG = East of Greenwich (degrees 
WOG = West of Greenwich (degrees) 
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At 45 degrees North and South parallel, the difference between the 
semi-major and semi-minor axis, as stated above, must be one-half of this 
value. It is 34. 5 m. 

At the poles, it is seen that the radial vector at the North Pole is longer 
than at the South Pole. The difference in height is 37. 8 m, which accounts for 
the "pear" shape of the earth, suspected already in 1959 by O'Keefe [9] . The 
gravity at the North Pole is higher by 13.11 mgal than at the South Pole. 

Around the equator belt, the gravity varies by as much as 31. 74 mgal (between 
the longitudes 72 degrees and 162 degrees east) . The mean equatorial gravity 
must be the same as for the triaxial ellipsoid, Model 4a; i.e. , 

g g = 978. 0320 gal . 

Model 5b ("Best" Geoid witho ut Sectorial Term) . For purposes of com- 
parison with models featuring circular equators, it is also of interest to deter- 
mine the shape of the "best" geoid, with all available zonal harmonics coeffi- 
cients, based on a circular equator. Model 5b is therefore obtained by solving 
equation (57) without the last (sectorial) term. Since the sectorial harmonic 
at the equator has already been determined in Figure 1, the desired circular- 
equator shape is obtained by subtracting it from the radius vectors of Model 5a. 


COMPARISON OF MODELS 


In accordance with the different assumptions used, the various geo- 
graphical theories to be compared here define geometrical bodies which differ 
from each other and from the geoid. Also, they lead to different gravity values 
for identical locations on the body. Primary attention is directed towards the 
deviations in geometries . 

Figure 6 shows the radius vector as a function of latitude of Models 1, 

2, and 3, referenced to the same body, i. e. , the "best" geoid without the 
sectorial term (Model 5b) . In all three cases, the bodies are not ellipsoids, 
but spheroids of revolution symmetrical around the equator plane, since they 
are based on superimposed spherical harmonics. It is seen that the widest 
deviation from the geoid is exhibited by Model 1, which takes into account 
solely the second order zonal harmonic, J 2 . At the North Pole, the difference 
is -32.6 m; at the South Pole, +5.20 m. Since the spheroid of Model 1, as 
shown on Table I, differs from the Model 1 ellipsoid (flattening 1/298.254) by 
as much as -18. 7 m at the North and South Pole, the deviation of the Model 1 


44 






ellipsoid from the reference geoid is less at the North Pole, amounting to -13. 9 
m, and higher (23. 9 m) at the South Pole. 

The spheroids of Models 2 and 3 are based on the old and the new set of 
coefficients J 2 , J3, J4, respectively, as given by equations (39) and (44). In 
general, they agree very well with the reference geoid, exhibiting deviations of 
less than 5 m. At the poles, the deviation of the Model 2 spheroid from the 
reference line as well as from the improved Model 3 spheroid is larger, pri- 
marily due to the changes in J 2 and J3. 

As in the case of Model 1, it is of interest to determine the deviation 
of the Model 2 spheroid from its associated ellipsoid, as given in Table II. The 
height of the spheroid based on the Standard values of J 2 , J 3 , J4 over the ellipsoid 
with consistent flattening 1/298.30 is shown in Figure 7. Also plotted is the 
height of the "best" geoid without the J 2j 2 -term (Model 5b) . The agreement 
between the two Model-2 bodies is good; the largest difference is at the South 
Pole, with 16.4 m. As is to be expected, the spheroid agrees better with the 
"best" geoid than the ellipsoid. The influence of particularly the odd zonal 
harmonics becomes manifest at the South Pole, where the Standard Ellipsoid, 
as used presently by NASA, deviates by as much as 27. 2 m from the Model 
5b-geoid. 

A still better approximation to the "best" geoid with circular equator 
can be given with the spheroid of Model 3, as depicted in Figure 8. Here, the 
spheroid height is referenced to the associated ellipsoid of flattening 1/298.222, 
derived from the set of improved J 2 , J3, and J 4 . The deviations of the spheroid 
from the ellipsoid are numerically very similar to those between the Model 2 
bodies, while the agreement of the ellipsoid with the circular-equator geoid has 
been improved considerably, except for the latitudes around 20N and 30S, where 
the gap is still almost 10 m. The agreement of the spheroid with the "best" 
geoid is generally better than 5 m. 

The height of the "best" geoid (Model 5a) has been computed from the 
data of Tables IV and V, and is presented in Figure 9 (greatly exaggerated) , 
referenced to the triaxial ellipsoid of Model 4a. Both bodies are based on 
expansions of spherical harmonics up to and including J 14 , with the latest 
available values for the coefficients . The sectional cut is along the prime 
meridian. The "pear" shape of the earth is clearly discernible. Since the 
mean equatorial radius is the same for both bodies (as well as the sectorial 
harmonic) , they share the same equator. At the poles, however, deviations 
of the geoid are observable, amounting to 18. 87 m at the North Pole, and to 
-18. 87 m at the South Pole. In the intermediate latitudes, the deviations are 
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North Pole 



Equator 

! 



South Polo 


FIGURE 8. HEIGHTS OF MODEL 3-SPHEROID AND "BEST" GEOID WITHOUT SECTORIAL TERM 
OVER ELLIPSOID OF FLATTENING 1/298.222 (MODEL 3) 









FIGURE 9. HEIGHT OF "BEST" GEO ID (MODEL 5) OVER TRIAXIAL 
ELLIPSOID OF FLATTENING 1/298. 184, ALONG PRIME MERIDIAN 




less, with -8. 5 m at 20N and 7. 7 m at 30S. The height differences are also 
plotted in Figure 10, on a rectilinear grid, for a meridional section at Greenwich 
longitude . 

According to the different harmonics coefficients used in the potential 
equation, the gravity of the analyzed models differs in each case. Using the 
gravitation of the "best" geoid (Model 5a) as tabulated in Table V, as common 
reference, the gravities of each model were computed for a meridional section 
along the prime meridian and plotted in Figure 11. As is to be expected, the 
gravity of the triaxial ellipsoid (Model 4a), listed in Table IV, agrees best 
(generally by better than ±3 mgal) with the "best" geoid. The other models, 
with gravity values determined from differing values and numbers of oblateness 
coefficients, show generally identical behavior, with significant deviations 
( up to 35 mgal) at the poles and equator, in accordance with the recent findings 
of a more pronounced influence of the odd zonal harmonics, particularly J 3 . 

For Model 1, the exclusive use of J 2 apparently compensates somewhat 
around the North Pole for the deviations in Models 2 and 3, which are obviously 
caused by the third order zonal harmonic, the only additional additive amplitude 
at these latitudes besides J 2 . At the lower latitudes, Model 3 — with the improved 
set of coefficients — agrees best of all three models with the reference, as is 
to be expected. However, the improvement is only very minor (<1 mgal) , 
leading to the observation that the effects of improving the oblateness coefficients 
(rather than increasing their number beyond J 4 ) are apparently revealed more 
significantly in the geometries than in the gravities of the bodies used. 


CONCLUSIONS 

1. A new earth figure model, or geoid, has been developed [2,3] and 
computed here, consistent with the most recent oblateness coefficients available. 

2. The currently used NASA Standard ellipsoid differs from this "best" 
geoid by -10.6 m at the North Pole, 9. 6 m at 20 degrees north parallel, -4. 6 m 
at 35 degrees south parallel, and as much as 27.2 m at the South Pole. If a 
spheroid based on a series expansion of spherical harmonics with the same 
oblateness coefficients as in the gravity equation were used instead, as provided 
by Krause's theory, these deviations could be reduced to 2. 4 m, 4. 9 m, 4. 0 m, 
and 10. 8 m, respectively. This surface formula would be more consistent with 
the gravity formulation. 
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South Pols 


FIGURE 10. HEIGHT OF "BEST" GEOID (MODEL 5) 
OVER TRIAXIAL ELLIPSOID OF FLATTENING 
1/298. 184, ALONG PRIME MERIDIAN 





FIGURE 11. DEVIATION OF GRAVITIES OF MODELS 1, 2, 3 AND 4 
FROM GRAVITY OF "BEST” GEOID (MODEL 5), 

ALONG PRIME MERIDIAN 



3. The currently used NASA Standard oblateness coefficients J 2 , J3, 

,J 4 result in gravity values, which differ from the gravity based on recent 
satellite measurements by as much as 35 mgal at the North Pole, -15. 7 mgal 
at the equator, and 30 mgal at the South Pole (along the Meridian) . Replacing 
the three J-coefficients by their updated values does not lead to significant 
improvements . 

4. Conversely, updating the oblateness coefficients J 2 , J3, J4 in the 
NASA Standard geometry model (and using the series expansion, as noted 
above) could lead to significant improvements in height deviations (Fig. 6) , 
amounting to about 100 percent at the North Pole, 50 percent at the South Pole, 
and 30-50 percent at most of the intermediate latitudes. 

5. For certain applications (e. g. , guidance of short duration first 
stage boosters) , use of geographical and gravitational models based only on 
J 2 may be permissible. Because of the shape of the J 2 -harmonic (Fig. 2), it 
is advisable to use the ellipsoid (flattening 1/298.254) , not the spheroid, as 
geometrical reference. Its height deviations from the "best" geoid (circular 
equator) are less than 10 m in the belt between 65 degrees north and 65 degrees 
south latitudes. The gravity of this simple model differs by maximally 30 mgal, 
at the South Pole (-18.5 mgal at equator, 16. 8 mgal at 90 degrees north) from 
the sophisticated gravity model. The deviation is zero at about 24 degrees 
north and 30 degrees south parallels intersecting the Prime. 

6. The latest available satellite measurements result in a "best" geoid 
which is pronouncedly "pear" shaped, with the north polar radius being 37. 8 m 
longer than the south polar radius vector. The equator is slightly elliptic, with 
a difference of 69. 0 m between the semi-major and semi -minor axes. The 
major axis of the equator ellipse is rotated West against the Prime Meridian. 

The rotation angle (longitude) is currently thought to be -18 degrees ±3 degrees. 
As shown by Krause, the ’’best" geoid can be described by a series expansion 

of spherical harmonics of the same degree as in the geopotential, the (amplitude) 
coefficients of which are directly related to the satellite -obtained oblateness 
coefficients . 
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7. As an approximation to the "best" geoid, an exact ellipsoid can be 
defined, as is customary in geodesy. The ellipsoid is a triaxial ellipsoid with 
the same equator ellipse as the geoid. The largest deviations of this body from 
the geoid are 18. 87 m at the poles. In a band bounded by the parallels 63 degrees 
north and 70 degrees south, the height differences amount to less than 10 m. 

If for the geopotential of the triaxial ellipsoid only even J's are used, the gravity 
deviates from the geopotential of the geoid generally by less than 7 mgal and 
between 70 degrees north and 70 degrees south by less than 5 mgal. 


George C. Marshall Space Flight Center 

National Aeronautics and Space Administration 

Marshall Space Flight Center, Alabama 35812, January 3, 1969 
933-50-07-01 
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APPENDIX A 


Computation of Legendre Polynomials 

The Legendre spherical functions P (x) are polynomials of the type 

n n(n-l) n-2 
_ X “ 2( 2n-i) X 

n(n-l) (n-2) (n-3) n-4 

+ 2 . 4(2n-i) ( 2n-3) X 

n = 0, 1, 2, . . . 

They are solutions of the differential equation 

( 1 - x 2 ) -^r%T - 2x ^ + n(n+l)y = 0. 

v 7 dx^ dx 

For the present analyses, the Legendre functions were computed from 
the recurrence relation 

(n+l)P n ^(x) - (2n + l)xP n (x) + nP^^fx) = 0, 

where 

x = sin0 
and starting with 

P 0 (si n<t>) = 1 
P^sin <p) = sin 0 . 

The most extensive tabulation of P n (cos0) , for n = 0( 1) 80 and 
c p = 0° ( 1° ) 180°, has been published in Reference 19. Since 

cos 0 = sin( 90° ± 0) , 

the spherical functions can also be taken from this reference. 


p n (x) = 


1 • 3 • 5. . . ( 2n-3) ( 2n-l) 


n: 
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APPENDIX B 


Alternate Expression for the Gravity Acceleration 

Krause, in References 2 and 3, gives a form of the gravity, 



oo 

C 0 " Yj (n + 1) C n P n (sin <p) + 3C 2> 2 cos 2 0 cos 2( A - A q ) , 
n=l 


obtained from equation (26) by several transformations [2]. The coefficients 

C are related to the oblateness coefficients J and can be determined from 
n n 

them by relations similar to the equations ( 16) through (22) for the coefficients 

A^, involving some of the latter. 


1 + 2(1 - A 0 ) - g cu e + — ( 3J 2 + co e ) 2 


1 


" g B 2 ( i^J 2 + 0> e 3B 2 ) 


Cl = 0 


J 2 - |a 2 " fw e “ £3 (3J 2 + £ e ) 2 + ^-B 2 (18J 2 + w e -3B 2 ) 


Ca = — 


J 4 - f A 4 + j|^-(3J 2 +w e ) 2 + ^B 2 (18J 2 + ^ e -3B 2 ) 


C = - 

K X 


2A 

K 


K + 1 


for k #0,2,4 


C 2, 2 = -J 2 ,2 


For the "best" geoid (Model 5a) with Kozai's set of J up to J 14 , the 
C -coefficients become n 
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1.001760646 


C 0 = 

C t = 0 

C 2 = - 0. 1172324264 x 10 -2 
C 3 = - 0.1272995717 x 10 -5 
C 4 = + 0.2636592115 x 10 -6 
C 5 = - 0. 1401281566 x 10" 6 
C 6 = + 0.4619357542 x 10“ 6 
C 7 = - 0.2500551096 x 10~ 6 
C 8 = - 0.2102749240 x 10 -6 
C 9 = -0.4245829102 x 10“ 7 
C 10 = - 0.4424481398 x 10“ 7 
C u = + 0.2520357773 x 10 -6 
C 12 = - 0.3025300610 x 10“ s 
C 13 = - 0.9786358495 x 10 -7 
C 14 = +0.1553740187 x 10 -6 

With these coefficients, equation (59) then becomes 6 , using 

g = 978.0320, 
e 

g = 979753.968 + 3439. 712P 2 ( sin<£) + 4. 980P 3 (sin <p) 

- 12. 893P 4 (sin0) + O.822P s (sin0) - 3. 163P e ( sirup) 

+ 1. 956P 7 (sin0) + 1 . 851P g ( sirup) + 0. 415P 9 ( sin0) 

+ 0. 476P lo (sin0) - 2. 958P U ( sirup) + 3. 846P 12 (sin0) 

+ 1. 34OP 13 (sin0) - 2. 279P 14 (sin<£) 

+ 15.873cos 2 0 cos2(X+18°) (mgal) . (60) 

6. Some of the amplitudes (coefficients) of the gravity equation in Reference 3 
have been found to be slightly incorrect. The values are corrected here. 
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With the Legendre polynomials obtained as shown in Appendix A, the 
gravity can be computed for the "best" geoid from equation ( 60) as an 
alternative to the gravity equation (26) . 
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